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ABSTRACT 

With  a  finite  volume  approach,  a  flux-form  semi-Lagrangian  (TFSL)  scheme  with  space-time  transformation  was  devel¬ 
oped  to  provide  stable  and  accurate  algorithm  in  solving  the  advection-diffusion  equation.  Different  from  the  existing  flux- 
form  semi-Lagrangian  schemes,  the  temporal  integration  of  the  flux  from  the  present  to  the  next  time  step  is  transformed  into 
a  spatial  integration  of  the  flux  at  the  side  of  a  grid  cell  (space)  for  the  present  time  step  using  the  characteristic-line  concept. 
The  TFSL  scheme  not  only  keeps  the  good  features  of  the  semi-Lagrangian  schemes  (no  Courant  number  limitation),  but  also 
has  higher  accuracy  (of  a  second  order  in  both  time  and  space).  The  capability  of  the  TFSL  scheme  is  demonstrated  by  the 
simulation  of  the  equatorial  Rossby-soliton  propagation.  Computational  stability  and  high  accuracy  makes  this  scheme  useful 
in  ocean  modeling,  computational  fluid  dynamics,  and  numerical  weather  prediction. 
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1.  INTRODUCTION 

From  a  physical  point  of  view,  advection  of  a  passive 
tracer  is  the  simple  transition  of  a  quantity  without  diffusion 
and  dispersion.  Numerical  approaches  in  atmospheric  and 
oceanic  modeling  inevitably  introduce  diffusion  (or  dissi¬ 
pation)  and  dispersion  into  the  approximate  solution.  The 
numerical  diffusion  and  dispersion  are  aliens  to  the  process 
that  is  being  modeled  (Chu  and  Fan  1998,  1999).  As  applied 
to  a  constituent  advection  problem,  these  numerical  artifacts 
manifest  themselves  as  nonphysical  mixing  by  numerical 
diffusion,  nonphysical  highs  and  lows  in  the  constituent  field 
caused  by  dispersion,  and  nonphysical  tracer  spectra  caused 
by  trapping  in  nonpropagating  small  spatial  scales  (Rood 
1987).  For  example,  the  commonly  used  upwind  scheme  is 
conditionally  stable  (with  the  Courant  number  being  much 
smaller  than  1)  and  some  artificial  viscosity  is  introduced. 
Hence,  less  the  numerical  diffusion  and  dispersion  errors 
equates  to  better  model  performance. 

Many  numerical  algorithms  have  been  proposed  to  re¬ 
duce  numerical  diffusion  and  dispersion  errors  and  to  keep 
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the  numerical  stability.  The  flux-form  semi-Lagrangian 
scheme  is  among  them.  Using  the  flux-form  semi-Lagran¬ 
gian  schemes,  artificial  viscosity  is  reduced  and  stability  is 
kept  without  the  limitation  of  a  Courant  number  (Casulli 
1990,  1999).  In  this  study,  we  use  a  finite  volume  approach 
to  develop  time-space  transformed  flux-form  semi-Lagran¬ 
gian  (TFSL)  scheme.  This  scheme  has  an  explicit  form  and 
much  less  diffusion  and  dispersion  errors. 

The  stability  and  accuracy  of  numerical  schemes  for 
ocean  models  are  usually  verified  using  the  propagation  of 
a  Rossby  soliton  on  an  equatorial  beta-plane.  In  principle, 
the  soliton  propagates  to  the  west  at  a  fixed  phase  speed, 
without  a  change  of  shape.  Since  the  uniform  propagation 
and  shape  preservation  of  the  soliton  are  achieved  through 
a  delicate  balance  between  linear  wave  dynamics  and  non¬ 
linearity.  In  other  words,  the  Rossby  soliton  is  non-diffusive 
and  non-dispersive  (Boyd  1980),  which  makes  it  a  perfect 
test  case  for  verification  of  numerical  schemes  in  ocean 
models  since  any  diffusion  and  dispersion  in  the  numeri¬ 
cal  solution  of  the  Rossby  soliton  are  computational  errors. 
Interested  readers  are  referred  to  the  website:  http://marine. 
rutgers.edu/po/index. php?model=test-problems. 
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To  show  the  benefit  of  using  the  TFSL  scheme,  we  first 
show  instability  and  large  diffusion  and  dispersion  errors  in 
numerical  solution  of  the  Rossby  soliton  using  the  existing 
schemes  such  as  the  flux-form  upwind,  flux-form  central, 
Lax-Wendroff,  and  flux-form  semi-Lagrangian  schemes. 
Then,  we  will  describe  the  procedures  of  the  TSFL  scheme 
development  and  verification.  The  rest  of  paper  is  organized 
as  follows.  Section  2  describes  the  equatorial  Rossby  soliton 
and  its  usefulness  for  an  ocean  model  verification.  Section 
3  shows  the  failure  of  the  three  existing  schemes  (upwind, 
central,  Lax-Wendroff,  semi-Largangian)  in  simulating  the 
equatorial  Rossby  soliton.  Section  4  introduces  the  TFSL 
scheme.  Section  5  derives  the  analytical  form  of  the  am¬ 
plification  factor  of  the  TFSL-scheme.  Section  6  shows  the 
capability  of  the  TFSL-scheme  in  simulating  the  equatorial 
Rossby  soliton.  Section  7  presents  our  conclusions. 


v(s,  y,  t)  =  2 y 
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and  the  variable  7)  (s,  t)  satisfies 
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f  =  1.5366,  /2  =  0.098765 
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which  is  the  Korteweg-de  Vries  (KDV)  equation  with  the 
exact  solution. 


2.  ROSSBY  SOLITON 

Let  Q  be  the  angular  frequency  of  earth’s  rotation  and 
R  be  the  earth  radius,  and  let  (x,  y)  be  the  spatial  coordinates 
with  unit  vectors  (i,  j)  and  t  be  the  time.  Consider  a  single 
layer  of  homogeneous  ocean  layer  with  depth  of  H.  Lamb’s 
parameter  is  defined  by 


4  Q2R2 
~  §H 


0) 


where  g  is  the  gravitational  acceleration.  The  length  and 
time  are  non-dimensionalized  by 


L  = 


Eui 

2  Q. 


(2) 


For  the  mean  ocean  depth  H  =  4  km,  the  earth  radius  R  = 
6370  km,  and  Q.  =  2ji  I  (86400  s),  the  length  and  time  scale 
are  L  =  543  km,  T  =  22.39  hr.  Let  (x,  y)  be  the  non-dimen¬ 
sional  Cartesian  coordinates,  ( u ,  v)  be  the  non-dimensional 
velocity  components  in  the  meridional  and  latitudinal  direc¬ 
tions,  and  <p  be  the  non-dimensional  surface  elevation.  After 
defining 


7j(s,  t)  =A  sech2[B(s  +  uB2t)] 

A  =  Q.T72B1,  B  =  0.394,  /J  =  0.395  (6) 


Substitution  of  the  exact  solution  Eq.  (6)  into  the  third  term 
in  the  lefthand  side  of  Eq.  (5)  leads  to 


d  7i 
dt 


=  s 


(7) 


S=f2 


d  3  7i 
8  s 3 


(8) 


where  S  is  treated  as  a  source/sink  term.  Evidently  Eq.  (7) 
has  the  analytical  solution  Eq.  (6).  Since  the  analytical  solu¬ 
tion  Eq.  (6)  exists,  the  Rossby  soliton  Eq.  (7)  is  a  perfect  test 
case  for  verifying  the  stability  and  accuracy  of  numerical 
schemes  since  the  diffusion  term  has  been  changed  into  the 
given  source/sink  term.  To  do  so,  the  fluid  is  assumed  to  oc¬ 
cupy  the  equatorial  region  surrounding  the  earth.  The  zonal 
direction  is  discretized  into  120  cells  (i.e.,  resolution  at  3° 
longitude).  The  increment  A,v  is  given  by 


A  s  = 


2ttR 
120  L 


0.256 


(9) 


s=x-  ct  (3) 

and  transforming  the  nonlinear  shallow  water  wave  equa¬ 
tions  into  a  frame  of  reference  moving  with  the  linear  wave, 
the  flow  variables  (m,  v,  0)  for  the  mode-1  wave  can  be  rep¬ 
resented  by  (Boyd  1980) 

(6y2  -  9)  /  y2\ 

m(s,}',0  =  — 4 - V(s,t)exp[-^-)  (4a) 


The  independent  variables  (,v,  t)  are  discretized  by  s,  =  s,  _ ,  + 
As,  t"  =  t"1  +  At,  i  =  1,  2,  ....;  n  =  1,  2, ....,  with  Af  the  time 
step.  The  dependent  variable  (>])  at  Ov, ,  t")  is  represented  by 
rfi  m  Vis,,  t"). 

3.  SEVERAL  EXISTING  SCHEMES 

Equation  (7)  can  be  discretized  using  the  flux-form  up¬ 
wind  scheme, 
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vT'  =  d"  +  )2-(-^')2]  +  2"  (10) 

the  flux-form  central  scheme, 

if*'  =  if  +  |^[(?7i'+i/2)2 -  W-w)2]  +  Q"  (11) 

the  Lax-Wendroff  scheme, 

V"*'  =  if  -  -  Vh)  +  -  2?7"  +  777-0  +  2"  (12) 

and  the  semi-Lagrangian  scheme, 

if*'  =  7)7  +  -  vZn)  +  2"  (13) 


so,  the  surface  elevation  0(s„  yy  t„)  is  plotted  with  contour 
values  of  2.13,  4.26,  6.4,  8.53,  10.66,  12.79,  14.93,  and 
17.06  cm.  All  the  numerical  schemes  greatly  distort  the 
Rossby  soliton  (Fig.  1).  When  the  ratio  of  the  root-mean 
square  error  versus  the  root-mean  of  the  analytical  solution 
is  greater  than  100%,  the  numerical  solution  is  considered 
divergence.  Figure  1  shows  that  the  numerical  solution  di¬ 
verges  at  27°45’W  using  the  flux-form  central  scheme,  at 
54°45’W  using  the  Lax-Wendroff  scheme,  and  at  30°W  af¬ 
ter  one  cycle  around  the  earth  using  the  flux-form  upwind 
scheme.  Comparing  Figs. lb  -  d  to  Fig.  la,  the  numerical 
solutions  are  totally  different  from  the  analytical  solution. 

4.  TFSL-SCHEME 

4.1  Semi-Lagrangian  Method 

Consider  the  advection  of  a  passive  scalar  0(x,  t)  by 
the  velocity  u(x,  t).  The  Eulerian  formulation  is  given  by 


with 


D(j) 

Dt 


d  0 
dt 


+  u  •  V0  =  S 


(16) 


2 A/2 

U 


[2sec/z2(T)  -  3 sec/z4(y)]|T 


(14a) 


where  x  is  the  position  vector,  D/Dt  denotes  the  material 
derivative,  while  the  Lagrangian  counterpart  is 


=  fv'lht 
As  ’ 


Y  =  B(s+pB2t) 


(14b) 


d$>„ 

dt  ~  'V 


dx„ 

dt 


=  u(x„,  t) 


(17) 


V’lta  =  j[(d'Um)2  +  (Vt*L)2],  Tftm  =  lf*vi-C,tW  (14c) 

In  order  to  compare  the  difference  between  numerical 
and  exact  solutions  (a  westward  propagating  Rossby  soli¬ 
ton),  the  zonal  equatorial  strip  is  assumed  to  be  infinitely 
long.  When  the  Rossby  soliton  travels  over  n  x  120  cells, 
it  goes  around  the  earth  once  n  times  (called  n  cycles).  The 
exact  solution  at  t  =  0  is  taken  as  the  initial  condition, 

7]  =  (s,  0)  =  A  sec  h2 {Bs)  (15) 

with  s  =  0  denoting  0°  longitude. 

Three  difference  Eqs.  (10)  -  (12)  are  solved  numerically 
from  the  initial  condition  Eq.  (15)  representing  the  upwind, 
central,  and  Lax-Wendroff  schemes  (Lax  and  Wendroff 
1960)  with  varying  At  at  each  time  step  for  a  given  Courant 

CA? 

number (C  =  0.75),  At  — - 77 - rr.  Selection ofC  =  0.75 

/|max(|  if  |) 

is  due  to  the  fact  that  the  proposed  TFSL  scheme  will  be  re¬ 
duced  to  the  Lax-Wendroff  scheme  for  C<0.5  [see  Eq.  (36)]. 
After  obtaining  the  numerical  solution,  77  (x„  f„),  substitut¬ 
ing  it  into  Eq.  (4c)  yields  0(s„  yy  f„).  The  accuracy  of  the 
schemes  can  be  verified  through  their  capability  in  predict¬ 
ing  the  westward  propagation  of  the  Rossby  soliton.  To  do 


where  the  subscript  ‘p’  shows  the  fluid  particle  in  Lagran¬ 
gian  sense.  Although  Eqs.  (16)  and  (17)  carry  the  same 
physical  information,  their  discretization  and  numerical 
implementation  is  different:  Eq.  (16)  is  discretized  on  an 
Eulerian  grid  with  a  finite  number  of  grid  points  and  then 
time-advanced,  while  Eq.  (17)  is  integrated  for  a  finite  num¬ 
ber  of  fluid  particles. 

Semi-Lagrangian  methods  combine  both  Eulerian  and 
Lagrangian  points  of  view;  the  scalar  field  is  discretized  on 
an  Eulerian  grid,  but  is  advanced  in  time  using  Eq.  (17).  The 
key  element  in  accomplishing  this  is  the  identification  of 
each  grid  point  x,  as  the  arrival  point,  for  instance,  at  f  +  At, 
of  a  particle  originating  from  x*  at  time  t.  The  algorithm  has 
three  steps:  (a)  The  particle  associated  with  each  grid  point 
x,  at  time  t  +  At  is  traced  back  to  its  location  x*  at  time  t, 

x*  =  x;-jT  u(r)dT  (18) 

(b)  The  scalar  value  at  (x* ,  t)  is  obtained  by  interpolating  the 
known  values  at  neighboring  grid  points, 

0(x*,f)  =  R[0fe),f]  (19) 
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where  P  is  any  interpolation  operator  and  (x,,)  denotes  the 
set  of  interpolation  points  associated  with  x*,  for  example, 
the  nodes  of  the  cell  containing  x* ;  (c)  Finally,  the  scalar  is 
updated, 

<$(x„  t  +  At)  =  <$(x*,  t)  +  Q,  (20) 

Thus,  the  main  issues  of  the  semi-Lagrangian  method  are 
the  backward  integration  in  step  (a)  and  the  interpolation 
in  step  (b). 

4.2  Flux  Form 

Equation  (16)  can  be  rewritten  in  the  flux  form  with 
inclusion  of  diffusion. 


-^y-  =  V-F  +  S,  F=-u  (p  +  KVcp  (21) 

where  k  is  the  diffusion  coefficient.  Let  the  dependent  vari¬ 
able  0(x,  t)  be  defined  on  the  space  0 ,  0  <  x  <  L, ,  0  <  y  < 
Ly.  0  <z<Lz- 

with  (L„  Ly,  L.)  the  lengths  in  (x,  y,  z)  directions.  Let 
Ax  =  ly  ,  Ay  =  jy  ,  Az  =  jf-  be  the  uniform  spatial  incre¬ 
ments  with  ( Nx  +l,Ny+  1,  Nz  +1)  the  grid  numbers.  Integrat¬ 
ing  Eq.  (21)  for  the  finite  volume,  A Q,jt  =  [x, .  m  <  x  <  xi+m, 

i  —  ,  Ax 

yj -in  —  y  —  Yj+ m,  Zk -m  —  z  —  Zk  +  mi,  Xi  +  m  —  xt+  ^  >»’±  1/2 

Ay  A" 

=  yj  ±  -j- ,  Zk±  112  =  zA  +  "2“,  from  tn  to  tn+  ,  =  tn  +  A  t, 

we  obtain  the  finite  difference  equation  of  the  flux-averaged 
transport, 


(a)  Analytic  Solution  with  Courant  Number=  0.75 


Fig.  1 .  Surface  elevation  (f>(s,  y,  t)  of  the  Rossby  solitons  obtained  from  an  (a)  exact  solution,  and  numerical  integration  with  C  =  0.75  using  the  (b) 
flux-form  upwind  scheme,  (c)  flux-form  central  scheme,  and  (d)  Lax-Wendroff  scheme.  Note  that  the  numerical  solution  diverges  at  30°W  after  one 
cycle  using  the  flux-form  upwind  scheme,  at  27°45’W  using  the  flux-form  central  scheme,  and  at  54°45’W  using  the  Lax-Wendroff  scheme. 
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j,  k~  0i,  j,  k 

At 


y7>(?t,  n  +  1)  T7>(n,  w  +  1)  S~i(n,n  + 1)  s~i(n,n  + 1) 

r  i  +  1/2,  j,k  -  r  i-  1/2,  y,  it  O'  i,  y  + 1/2,  it  -  O  i,  j  -  1/2,  it 

Ax  +  Ay 

7/(«,n  +  l)  7T(n,«  +  l) 

ii  y,  it +  1/2  "  -O  /,/,*- 1/2  A 

+ - £ - +s"" 


(22) 


In  the  existing  flux-form  semi-Lagrangian  schemes,  the 
temporally  integrated  flux  [similar  for  F  (+1/2  ]  is 

given  by  the  mean  value  at  the  two  time  steps  t„  and  t„  +  , 
(e.g.,  Casulli  1990;  Tanguay  et  al.  1990), 


where  ( F ,  G,  FT)  are  components  of  the  vector  F,  and 


77>(/i,  n  +  1) 

r  i  - 1/2 


=  \(Flm  +  FVm) 


(28a) 


=)(„.„+',  =  1  j  Fdf 


Using  the  characteristic-line  concept,  the  flux  at  time  step 
t„  +  j  and  location  x,  1/2  can  be  transformed  into  the  flux  at 
time  step  t„  and  location  x,-.  1/2.c  (Fig.  2), 


represents  the  temporal  average  (from  t„  to  t„  +  i).  The  tilde 

represents  the  volume  average  over  Qijt ,  F"*m  =  F".m-c  (28b) 


uk  =  fahyh?  fff  &  dxdydz  (24a) 

Qjk 

The  hat  represents  the  combined  volume  (£2tfi)  and  temporal 
average  (from  t„  to  t„  + 1), 

/„  + 1 

^  =  AtAxAyA?  f  Iff  S  dxdydzdt  (24b> 

'  r,  Ogt 


For  the  finite  volume  AQijt ,  the  flux  at  x  =  x,  _  ln  and  t  =  t„  is 
calculated  by 


Substitution  of  Eq.  (28b)  into  Eq.  (28a)  gives 

Ft^  =  ^{F-.m  +  Flm-c)  (28c) 

with  the  mean  flux  Ftm"  determined  at  the  current  time 
step  t„ .  Here 


_  u At 

Ax 


(29) 


is  the  Courant  number. 

Here,  we  propose  a  new  method  to  compute  the  tem¬ 
porally  averaged  flux  Ft m"  with  the  transformation  into 
spatial  averaged  flux. 


To  solve  Eq.  (22)  numerically,  we  need  to  compute  the  tem¬ 
porally  integrated  fluxes,  Ftimj.k,  FtvVlk,  Gu'/lm.k,  G 
II  /.  j,  k  + 1/2,  H  Ilk- 1/2.  If  these  fluxes  are  computed  using  the 
semi-Lagrangian  method,  it  is  called  the  flux-form  semi- 
Lagrangian  scheme  (Casulli  1990,  1999;  Lin  and  Rood 
1996). 


4.3  Transformation  of  Temporal  Integration  into  Spa¬ 
tial  Mean 

For  simplicity  and  no  loss  of  generality,  we  consider 
one  dimensional  problem  of  Eq.  (22)  without  source/sink 
term  (i.e.,  Sijt  =  0), 


i  rt,,+At  i  rxi.m 

FT-wl)  -  j[  F(Xi.m,  t ) dt  -  c'AtJ  'JF{x,  Q dx  (30) 


Substitution  of  Eq.  (29)  into  Eq.  (30)  leads  to 


Ff-t 


FI.  i/2  +  FI 


1/2  '  1  i-m-C 


2  PCf  2 

8m{Flm  +  FU)  £&(*7-*  +  FT.t.i) 

2  +  2 

ifC>^ 

(31) 


$*' -  $  Ftim"  -  Ftm" 

At  Ax 


(26) 


where 


From  the  semi-Lagrangain  consideration,  we  have 


C  1 
C  2 


,  <5„ 2  =  2^,  8t  =  ^r,  8m  =  1  -  5„ 


C’ 


(32) 


(x,*,  t„)  =  0  (x„  t„)  + 


[FJ! 


-  Ftml)]  At 

Ax 


The  bracket  [  ]  represents  the  round-off  integer.  Similarly, 
the  temporally  averaged  flux  at  the  right  boundary  (x  =  xi+  1C) 
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j-t  (n,  n  + 1) 

r  i  +  1/2 


F  i  + 1/2  +  F  i  + 1/2 .  c  -r  ^  1 

- 2 -  lf  C  —  2 

111  - 1 

(F"+1/2  +  E”)  ('Fi' *+l  +  F'-^ 


+  a-(^,-+^w.c)  jfoi 


(33) 


The  temporally  averaged  fluxes  Ft  "a"  and  F  S+1/2 0  (from  t„ 
to  t„  +  At)  are  transformed  into  the  spatially  averaged  fluxes 
over  multiple  grids  at  time  step  tn  with  weights  of  dm,  <5i, 
dm.  If  the  characteristic  line  at  t„  is  beyond  the  boundary,  the 
boundary  condition  can  be  used  to  calculate  Ftmn  (Fig.  3), 


l(F1  +  Fy2)  +  (c-l)(F"b  +  F") 

F& f  *>  =  - 2C  -  (34) 

where  F"b  is  the  boundary  value  between  F"  and  F'[ + 1 ,  and 
is  interpolated  by 


F(x,tn+At) 


n+1 


F(x-cAt,tn) 
t 

i-1  i— 1/2  i  i+1/2  i+1 

Fig.  2.  Temporally  varying  flux  at  the  boundary  x,_  y2  from  t„  to  t„  + 
At  is  transformed  into  spatially  varying  flux  at  f„  from  xt .  y2  -  CAt  to 
Xj  _  1/2  using  the  characteristic-line  concept. 


Fig.  3.  Same  as  Fig.  2  except  at  the  left  boundary  of  the  integration 
domain. 


n  =  (l-w)F'*1+i:F"  (35) 

Substitution  of  Eqs.  (31)  and  (33)  into  the  difference 
Eq.  (26)  leads  to 


Different  schemes  have  different  algorithms  to  compute  the 
temporally  averaged  fluxes  Ft  in  and  Ft  in'*  (from  tn  to  t„ 
+  At).  The  TFSL  scheme  has  second  order  accuracy  in  both 
time  and  space. 


0?  +  l  =  0"  + 


-  y  ( 0"+ 1  -  <t>l  1 )  +  y-(  0"+ 1  -  20"  +  0".  1),  C  <  j 
-j(<!>ui  -  07- 1)  -  y(0?  +  07-i  -  07 r,  -  07™ -1) 

-  (d  -  y- )(07™  -  07™-i)  -  y-(07™-i  -  07™ -2) 
+  -jjr(0"+,  -20.  +07i),  c  >  y 

(36) 


which  is  called  the  Transformed  Flux-formed  Semi-Lagran- 
gian  (TFSL)  scheme  for  the  advection-diffusion  Eq.  (21). 
Here,  D  =  C  -  m  -  1/2.  The  major  difference  between  the 
existing  flux-form  semi-Lagrangian  scheme  and  the  TFSL 
scheme  comes  from  the  different  calculation  of  the  tempo¬ 
rally  averaged  flux  Ft T/2+1):  Eq.  (28c)  for  the  existing  flux- 
form  semi-Lagrangian  scheme  and  Eq.  (31)  for  the  TFSL 
scheme. 

For  C  <  2,  the  TFSL  scheme  is  the  same  as  the  Lax- 
Wendroff  scheme.  Compared  to  the  central  difference 
(CED),  the  TFSL-scheme  has  an  extra  positive  term, 

TFSL-CED  =  -y  (0"+1  -  20"  +  070  (37) 


for  C  <  1/2.  This  term  can  be  regarded  as  the  numerical 
(positive)  diffusion  which  leads  to  computational  stability. 


5.  STABILITY  OF  THE  TFSL  SCHEME 

The  stability  of  numerical  schemes  is  an  important 
issue  in  solving  the  advection  Eq.  (16).  In  section  3,  we 
showed  the  instability  of  the  existing  schemes  (upwind, 
central,  and  Lax-Wendroff).  To  determine  the  stability  of 
the  TFSL  scheme  Eq.  (36),  the  Fourier  series  expansion  is 
used.  Decay  or  growth  of  an  amplification  factor  indicates 
whether  or  not  the  numerical  algorithm  is  stable  (von  Neu¬ 
mann  and  Richtmyer  1950).  Assuming  that  at  any  time  step 
tn,  the  compute  solution  0"  is  the  sum  of  the  exact  solution 
0"  l“°  and  error  e" , 

0?  =  0?W  +  £?  (38) 

and  substituting  Eq.  (38)  into  Eq.  (36),  we  obtain 

-y(e"+i  -  e7i)  +  y-(<s7-i  -  2e"  +  e7i),  C  57  2" 

-  e7i)  -  +  eh  -  e".m  -  e7™-i) 

S"  + 1  =  £”  + 

-  (d  -  y  )(e7 ™  -  £7™-i)  -  y-(e7„-i  -  e7™-2) 
+  -  2s”  +  Si- 1 ),  C  >  y 

(39) 
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The  finite  mesh  function,  £",  can  be  decomposed  into  a  Fou-  which  implies  that  the  TFSL-scheme  Eq.  (36)  is  stable  for 
ner  series,  all  the  C  values  (without  Courant  number  restriction). 


Nx 

e"  =  2  a"exp  (/;'#),  9  =  jiz/Nx  (40) 

j=-N, 

with  I  =  V- 1  ,  ( a ",  0)  being  the  amplitude  and  phase  an¬ 
gle  of  the  jth  harmonic.  Substituting  Eq.  (40)  into  Eq.  (39) 
yields 

a"*'  =  g(9,  C)a"  (41) 

where 


6.  SIMULATING  THE  ROSSBY  SOLITON  USING 
THE  TFSL  SCHEME 

The  TFSL-scheme  Eq.  (36)  is  only  for  a  spatially  vari¬ 
ant  and  temporally  invariant  u.  When  u  [or  -j\ rj  in  Eq.  (7)] 
at  x, .  ]/2  varies  with  time  from  t„  to  tn+1,  concept  of  variant 
characteristic  lines  can  be  used  to  determine  w(x,_  !/2,  t)  with 
sub  time-steps  (dtm,  dtx ...,  dtm  )  (between  t„  and  tn  +  i)  from 
u(x,  t„ )  at  grid  points  (x,.  i,...,  x,-  m,  x,»),  and  for  u  >  0  the 
time  from  the  left  neighboring  grid  x,  [k  +  ,  |  to  x, .  k  is  given 
by  (Fig.  5), 


g(d,C) 


1  -  C2(l  -  cos  tf)  -  7C  sin  ft  C<j 
^(1  -  cos  9)  +  [\  -  D  +  R^Jcos(mff) 

+  [j  +  D-D1)cos[(m-  1)9] 
+  44  cos  [( m  -  2)6 ] 


■  I]  [  y  -  D  +  ^  j  sin  (mO) 

+  {\  +  D-D2)sm[(m-l)9] 

+  ^sin[(w-2)ft]},  C>j 


(42) 


is  called  the  amplification  factor,  whose  magnitude  is  given 
by 


5,  _  7  dx  _  Ax  f  dz  _  A t  ln(l  - 

k~  J  u(x,t„)~ulkJ  1  -  Cl .tZ~~G.it  £ 

x,-l*  +  i]  0 

=  ~T7(1+%L  +  %L+-)’  ^  =  0.5, 1,2,..., m  (44) 


where 


A  .1  —  Mi.(t  +  1)  ^  It ,  :  At  /  A  C  \ 

Ax  =  X;.t  -  Xj_(t  +  1),  C.i  =  - (<7T - ,  C,.k  =  (45) 

The  parameter  C,.*is  the  Courant  number  for  sub  time  steps. 
A  formula  similar  to  Eq.  (44)  can  be  obtained  for  m  <  0  (us¬ 
ing  the  right  neighboring  grid).  The  temporally  averaged 
fluxes  from  f„  to  t„  +  Af  can  be  calculated  by  {taking  Ff:"a') 
[see  Eq.  (31)]  as  the  example} 


Is  (0,0  P  = 


[l-C2(l-costf)]2+C2sin2ft  C< 


(1  -  cos 


n  D 
2  ~D+  2 


+(^  +  z)-d2)2  +  ^ 

+  4(1  -  cos  #)  j(4  -  D  +  Rj-jcos (m6) 
4  +  D  -  D~  }cos  [(m  - l)ft] 


D2 


cos  [(hi 


-2)0]} 


+  (l-2D  +  2D1)[j  +  D-D2)cos9 


+  D 


ih 


D  +  RA  cos  (2(9),  C>4 


t7>(h,  n  +  1) 

r  i  - 1/2  = 


Ar 


o  Fi_vl  + Fj.i  £  F.i  +  F.i 

Ot  1/2  o  +  O^i  ^  +  ... 


2 

F,m  +  F, 


(46) 


Equation  (7)  for  the  Rossby  soliton  is  discretized  using 
the  flux  form. 


•Y)n  + 1  W'  T-i  (n,  n  +  1)  t-t(h,  n  +  1) 

//i  “  //;  r  (  +  1/2  -  r  i-m 


At 


A.? 


-  +  S, 


(47) 


where  .S',  is  the  temporally-spatially  averaged  source  term 


tn  +  1  S/  +  1/2 

Si  =  /  /  S(s,  t)dsdt 


(48) 


(43) 

The  TFSL-scheme  is  computationally  stable  if  |  g(9,  C)  <  1 
and  computationally  unstable  if  \g(9,C)  |  >  1.  Figure  4 
shows  that  g  (9,  C)  \  <  1  for  all  6  and  C  (larger  than  20), 


with  S(s,  t)  given  by  Eq.  (8).  The  difference,  Eq.  (47),  is 
solved  numerically  from  the  initial  condition  Eq.  (15)  using 
the  TFSL-scheme.  To  compare  with  the  existing  flux-form 
semi-Lagrangian  scheme,  the  Courant  number  is  set  to  1.5. 
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CFL  number 

Fig.  4.  Dependence  of  the  amplification  factor  g  (0,  C)  |  °f  the  TSFL  scheme  on  6  and  C. 


Wl 


<l 


i-m-1  i  i-m  i-2  i-1  i— 1/2 

Fig.  5.  Same  as  Fig.  2  except  for  temporally  varying  u. 


After  the  numerical  solution  i](xn  t„ )  is  obtained,  substitut¬ 
ing  it  into  Eq.  (4c)  yields  ([)(s„  y„  /„)  as  shown  in  Fig.  6. 
Note  that  the  flux-form  semi-Lagrangian  scheme  is  highly 
distorted  (Fig.  6)  with  the  numerical  solution  diverging  at 
40°30’W,  which  may  be  caused  by  the  error  accumulation. 
However,  the  TFSF-scheme  is  quite  stable  and  accurate. 
After  propagating  westward  around  the  earth  the  numeri¬ 
cal  Rossby  soliton  (using  the  TFSF  scheme)  appears  to  be 
almost  non-diffusive  and  non-dispersive. 

To  show  the  quality  of  the  TFSF-scheme,  the  differ¬ 
ence  Eq.  (47)  is  integrated  for  C-  1.5  for  a  long  time  period 
corresponding  to  the  Rossby  soliton  propagates  westward 
around  the  earth  5  times.  The  solution  0(s„  yn  t„ )  is  stable 
all  the  time  (Fig.  7).  The  relative  root-mean-square  error 
(rrmse), 


rrmse  (f)  = 


EE[r,(s!,yJ,t)-ri(si,fct)]2 


■'=i  i 


max  |  <piex)(s i,  y Jt  t) 


(49) 


is  calculated  to  illustrate  the  accuracy  of  the  TFSL  scheme. 
Table  1  shows  RRMSE  at  the  end  of  first  five  cycles  around 
the  earth.  The  error  varies  from  2.66%  for  the  first  cycle  to 
3.53%  for  the  fifth  cycle. 


7.  CONCLUSIONS 

(1)  This  study  shows  that  the  TFSL  scheme  is  a  promising 
stable  and  accurate  method  for  solving  the  advection- 
diffusion  equation.  The  Fourier  analysis  shows  that  the 
TFSL  scheme  has  second-order  accuracy  in  time  and 
space.  This  scheme  retains  the  good  features  of  semi- 
Lagrangian  schemes  (no  Courant  number  limitation) 
with  higher  accuracy.  Computational  stability  and  high¬ 
er  accuracy  than  the  widely  used  schemes  (central,  up¬ 
wind,  Lax-Wendroff,  semi-Lagragian)  makes  this  tech¬ 
nique  useful  in  ocean  modeling,  computational  fluid 
dynamics,  and  numerical  weather  prediction. 

(2)  Several  major  features  distinguish  the  TFSL  scheme 
from  existing  schemes,  both  Eulerian  and  semi-Lagran¬ 
gian.  First,  the  flux  (F)  at  the  side  of  each  grid  cell  is 
computed  not  from  a  single  time  step  (present  or  next) 
but  from  a  temporal  integration  from  the  present  time 
step  to  the  next  time  step.  Second,  this  temporal  integra¬ 
tion  is  transformed  into  a  spatial  integration  at  the  pres¬ 
ent  time  step  using  the  characteristic  line  method. 

(3)  The  equatorial  Rossby  soliton  is  used  to  test  the  capa¬ 
bility  of  the  TFSL  scheme  since  it  has  exact  solution. 
The  equation  is  solved  numerically  from  the  soliton  ini¬ 
tially  located  at  the  equator  and  0°  longitude  with  an 
overall  Courant  number  of  0.75.  The  upwind,  central, 
and  Lax-Wendroff  schemes  greatly  distort  the  Rossby 
soliton  and  diverge  as  it  propagates.  However,  the  TFSL 
scheme  does  not  distort  the  Rossby  soliton  and  converg¬ 
es  as  it  propagates  many  cycles  around  the  earth.  With 
an  overall  Courant  number  of  1.5,  the  numerical  Rossby 
soliton  can  still  propagate  many  cycles  around  the  earth 
using  the  TFSL  scheme,  but  diverges  at  40°30’ W  using 
the  flux-form  semi-Lagrangian  scheme. 

(4)  Application  of  the  TFSL  scheme  to  the  atmospheric  and 
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(a)  Analytic  Solution  with  Courant  Number=  1 .5 


Longitude 

Fig.  6.  Surface  elevation  (p(s,  y,  t)  of  the  Rossby  soliton  obtained  from  (a)  exact  solution,  and  numerical  integration  with  C  =  1.5  using  (b)  TFSL- 
scheme,  and  (c)  flux-form  semi-Lagrangian  scheme.  The  solutions  (p(s ,  y,  t)  are  plotted  at  four  time  instances  for  the  Rossby  soliton  (exact  solu¬ 
tion)  westward  propagating  90°,  180°,  270°,  and  360°  (return  to  the  initial  location). 


Cycle:  4  RRMSE:  2.37% 


Lontinude 


Cycle:  5  RRMSE:  2.41% 


-20  0  20  40 

Lontinude 


Fig.  7.  Surface  elevation  (p(s ,  y,  t )  of  the  Rossby  soliton  after  1-5  cycles  around  the  earth  obtained  from  numerical  integration  with  C  =  1.5  using 
the  TFSL  scheme. 
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Table  1.  RRMSE  of  the  surface  elevation  predicted  using  the  TFSL- 
scheme  after  the  first  five  cycles  around  the  earth. 


Cycle 

1 

2 

3 

4 

5 

RRMSE  (%) 

2.66 

2.86 

3.00 

3.22 

3.53 

oceanic  models  needs  more  research.  This  is  because 
that  the  highly  accurate  treatment  of  the  source  term  en¬ 
abled  good  performance  of  the  TFSL  scheme.  However, 
this  is  a  very  special  case  because  the  analytical  solution 
of  S  (source  term)  is  known  for  this  system.  In  ocean 
modeling,  computational  fluid  dynamics,  or  numerical 
weather  prediction,  source  term  analytical  solutions  are 
usually  unknown  and  thus  several  iteration  processes 
will  be  required  for  the  departure/arrival  point  estima¬ 
tion  as  well  as  the  source  term  estimation  (for  spatial 
and  temporal  averaging),  which  causes  the  loss  of  ef¬ 
ficiency  and  accuracy. 

(5)  The  TFSL  scheme  was  developed  on  the  basis  of  a  fi¬ 
nite  volume  approach.  It  is  relatively  easy  to  extend 
one-dimensional  space -time  transformation  Eqs.  (28b) 
and  (28c)  into  three-dimensional  transformations.  The 
space  integration  of  the  flux  is  conducted  over  the  two- 
dimensional  surface  of  the  finite  volume,  and  the  time 
integration  is  for  that  volume  (see  section  4.2).  This  will 
be  reported  in  a  separate  paper  in  the  near  future. 
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